This training manual provides a comprehensive guide to handling geospatial data in R, covering both vector and raster data formats, coordinate reference systems, and various mapping techniques.
By the end of this tutorial, you will be able to:
sf
packageterra packageggplot2 and
tmapleaflet and
mapviewFirst, let’s install and load all necessary packages:
# Install packages if not already installed
# install.packages(c("sf", "terra", "raster", "ggplot2", "tmap",
# "leaflet", "mapview", "dplyr", "viridis",
# "rnaturalearth", "rnaturalearthdata", "lwgeom"))# Load required libraries
library(sf) # Simple Features for vector data
library(terra) # Spatial data analysis (raster)
library(raster) # Legacy raster package
library(ggplot2) # Data visualization
library(tmap) # Thematic mapping
library(leaflet) # Interactive maps
library(mapview) # Quick interactive viewing
library(dplyr) # Data manipulation
library(viridis) # Color palettes
library(rnaturalearth) # World map data
library(rnaturalearthdata)
# Set tmap mode to plot (static) initially
tmap_mode("plot")Vector data represents geographic features as points, lines, or polygons with associated attributes.
The sf package is the modern standard for handling
vector geospatial data in R. It integrates seamlessly with the tidyverse
and provides a consistent interface.
# Create simple point features manually
points_data <- data.frame(
name = c("London", "Paris", "Berlin", "Madrid", "Rome"),
population = c(8982000, 2161000, 3645000, 3223000, 2873000),
lon = c(-0.1278, 2.3522, 13.4050, -3.7038, 12.4964),
lat = c(51.5074, 48.8566, 52.5200, 40.4168, 41.9028)
)
# Convert to sf object
cities_sf <- st_as_sf(points_data,
coords = c("lon", "lat"),
crs = 4326) # WGS84 coordinate system
# Display structure
print(cities_sf)## Simple feature collection with 5 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -3.7038 ymin: 40.4168 xmax: 13.405 ymax: 52.52
## Geodetic CRS: WGS 84
## name population geometry
## 1 London 8982000 POINT (-0.1278 51.5074)
## 2 Paris 2161000 POINT (2.3522 48.8566)
## 3 Berlin 3645000 POINT (13.405 52.52)
## 4 Madrid 3223000 POINT (-3.7038 40.4168)
## 5 Rome 2873000 POINT (12.4964 41.9028)
## Geometry set for 5 features
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -3.7038 ymin: 40.4168 xmax: 13.405 ymax: 52.52
## Geodetic CRS: WGS 84
## [1] POINT POINT POINT POINT POINT
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE
## xmin ymin xmax ymax
## -3.7038 40.4168 13.4050 52.5200
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
## [1] 5
## name population geometry
## Length:5 Min. :2161000 POINT :5
## Class :character 1st Qu.:2873000 epsg:4326 :0
## Mode :character Median :3223000 +proj=long...:0
## Mean :4176800
## 3rd Qu.:3645000
## Max. :8982000
# Get world country boundaries from Natural Earth
world <- ne_countries(scale = "medium", returnclass = "sf")
# Display first few rows
head(world)## Simple feature collection with 6 features and 168 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.36621 ymin: -22.40205 xmax: 109.4449 ymax: 41.9062
## Geodetic CRS: WGS 84
## featurecla scalerank labelrank sovereignt sov_a3 adm0_dif level
## 1 Admin-0 country 1 3 Zimbabwe ZWE 0 2
## 2 Admin-0 country 1 3 Zambia ZMB 0 2
## 3 Admin-0 country 1 3 Yemen YEM 0 2
## 4 Admin-0 country 3 2 Vietnam VNM 0 2
## 5 Admin-0 country 5 3 Venezuela VEN 0 2
## 6 Admin-0 country 6 6 Vatican VAT 0 2
## type tlc admin adm0_a3 geou_dif geounit gu_a3 su_dif
## 1 Sovereign country 1 Zimbabwe ZWE 0 Zimbabwe ZWE 0
## 2 Sovereign country 1 Zambia ZMB 0 Zambia ZMB 0
## 3 Sovereign country 1 Yemen YEM 0 Yemen YEM 0
## 4 Sovereign country 1 Vietnam VNM 0 Vietnam VNM 0
## 5 Sovereign country 1 Venezuela VEN 0 Venezuela VEN 0
## 6 Sovereign country 1 Vatican VAT 0 Vatican VAT 0
## subunit su_a3 brk_diff name name_long brk_a3 brk_name brk_group
## 1 Zimbabwe ZWE 0 Zimbabwe Zimbabwe ZWE Zimbabwe <NA>
## 2 Zambia ZMB 0 Zambia Zambia ZMB Zambia <NA>
## 3 Yemen YEM 0 Yemen Yemen YEM Yemen <NA>
## 4 Vietnam VNM 0 Vietnam Vietnam VNM Vietnam <NA>
## 5 Venezuela VEN 0 Venezuela Venezuela VEN Venezuela <NA>
## 6 Vatican VAT 0 Vatican Vatican VAT Vatican <NA>
## abbrev postal formal_en
## 1 Zimb. ZW Republic of Zimbabwe
## 2 Zambia ZM Republic of Zambia
## 3 Yem. YE Republic of Yemen
## 4 Viet. VN Socialist Republic of Vietnam
## 5 Ven. VE Bolivarian Republic of Venezuela
## 6 Vat. V State of the Vatican City
## formal_fr name_ciawf note_adm0 note_brk
## 1 <NA> Zimbabwe <NA> <NA>
## 2 <NA> Zambia <NA> <NA>
## 3 <NA> Yemen <NA> <NA>
## 4 <NA> Vietnam <NA> <NA>
## 5 República Bolivariana de Venezuela Venezuela <NA> <NA>
## 6 <NA> Holy See (Vatican City) <NA> <NA>
## name_sort name_alt mapcolor7 mapcolor8 mapcolor9 mapcolor13 pop_est
## 1 Zimbabwe <NA> 1 5 3 9 14645468
## 2 Zambia <NA> 5 8 5 13 17861030
## 3 Yemen, Rep. <NA> 5 3 3 11 29161922
## 4 Vietnam <NA> 5 6 5 4 96462106
## 5 Venezuela, RB <NA> 1 3 1 4 28515829
## 6 Vatican (Holy See) Holy See 1 3 4 2 825
## pop_rank pop_year gdp_md gdp_year economy
## 1 14 2019 21440 2019 5. Emerging region: G20
## 2 14 2019 23309 2019 7. Least developed region
## 3 15 2019 22581 2019 7. Least developed region
## 4 16 2019 261921 2019 5. Emerging region: G20
## 5 15 2019 482359 2014 5. Emerging region: G20
## 6 2 2019 -99 2019 2. Developed region: nonG7
## income_grp fips_10 iso_a2 iso_a2_eh iso_a3 iso_a3_eh iso_n3
## 1 5. Low income ZI ZW ZW ZWE ZWE 716
## 2 4. Lower middle income ZA ZM ZM ZMB ZMB 894
## 3 4. Lower middle income YM YE YE YEM YEM 887
## 4 4. Lower middle income VM VN VN VNM VNM 704
## 5 3. Upper middle income VE VE VE VEN VEN 862
## 6 2. High income: nonOECD VT VA VA VAT VAT 336
## iso_n3_eh un_a3 wb_a2 wb_a3 woe_id woe_id_eh woe_note
## 1 716 716 ZW ZWE 23425004 23425004 Exact WOE match as country
## 2 894 894 ZM ZMB 23425003 23425003 Exact WOE match as country
## 3 887 887 RY YEM 23425002 23425002 Exact WOE match as country
## 4 704 704 VN VNM 23424984 23424984 Exact WOE match as country
## 5 862 862 VE VEN 23424982 23424982 Exact WOE match as country
## 6 336 336 -99 -99 23424986 23424986 Exact WOE match as country
## adm0_iso adm0_diff adm0_tlc adm0_a3_us adm0_a3_fr adm0_a3_ru adm0_a3_es
## 1 ZWE <NA> ZWE ZWE ZWE ZWE ZWE
## 2 ZMB <NA> ZMB ZMB ZMB ZMB ZMB
## 3 YEM <NA> YEM YEM YEM YEM YEM
## 4 VNM <NA> VNM VNM VNM VNM VNM
## 5 VEN <NA> VEN VEN VEN VEN VEN
## 6 VAT <NA> VAT VAT VAT VAT VAT
## adm0_a3_cn adm0_a3_tw adm0_a3_in adm0_a3_np adm0_a3_pk adm0_a3_de adm0_a3_gb
## 1 ZWE ZWE ZWE ZWE ZWE ZWE ZWE
## 2 ZMB ZMB ZMB ZMB ZMB ZMB ZMB
## 3 YEM YEM YEM YEM YEM YEM YEM
## 4 VNM VNM VNM VNM VNM VNM VNM
## 5 VEN VEN VEN VEN VEN VEN VEN
## 6 VAT VAT VAT VAT VAT VAT VAT
## adm0_a3_br adm0_a3_il adm0_a3_ps adm0_a3_sa adm0_a3_eg adm0_a3_ma adm0_a3_pt
## 1 ZWE ZWE ZWE ZWE ZWE ZWE ZWE
## 2 ZMB ZMB ZMB ZMB ZMB ZMB ZMB
## 3 YEM YEM YEM YEM YEM YEM YEM
## 4 VNM VNM VNM VNM VNM VNM VNM
## 5 VEN VEN VEN VEN VEN VEN VEN
## 6 VAT VAT VAT VAT VAT VAT VAT
## adm0_a3_ar adm0_a3_jp adm0_a3_ko adm0_a3_vn adm0_a3_tr adm0_a3_id adm0_a3_pl
## 1 ZWE ZWE ZWE ZWE ZWE ZWE ZWE
## 2 ZMB ZMB ZMB ZMB ZMB ZMB ZMB
## 3 YEM YEM YEM YEM YEM YEM YEM
## 4 VNM VNM VNM VNM VNM VNM VNM
## 5 VEN VEN VEN VEN VEN VEN VEN
## 6 VAT VAT VAT VAT VAT VAT VAT
## adm0_a3_gr adm0_a3_it adm0_a3_nl adm0_a3_se adm0_a3_bd adm0_a3_ua adm0_a3_un
## 1 ZWE ZWE ZWE ZWE ZWE ZWE -99
## 2 ZMB ZMB ZMB ZMB ZMB ZMB -99
## 3 YEM YEM YEM YEM YEM YEM -99
## 4 VNM VNM VNM VNM VNM VNM -99
## 5 VEN VEN VEN VEN VEN VEN -99
## 6 VAT VAT VAT VAT VAT VAT -99
## adm0_a3_wb continent region_un subregion
## 1 -99 Africa Africa Eastern Africa
## 2 -99 Africa Africa Eastern Africa
## 3 -99 Asia Asia Western Asia
## 4 -99 Asia Asia South-Eastern Asia
## 5 -99 South America Americas South America
## 6 -99 Europe Europe Southern Europe
## region_wb name_len long_len abbrev_len tiny homepart
## 1 Sub-Saharan Africa 8 8 5 -99 1
## 2 Sub-Saharan Africa 6 6 6 -99 1
## 3 Middle East & North Africa 5 5 4 -99 1
## 4 East Asia & Pacific 7 7 5 2 1
## 5 Latin America & Caribbean 9 9 4 -99 1
## 6 Europe & Central Asia 7 7 4 4 1
## min_zoom min_label max_label label_x label_y ne_id wikidataid
## 1 0 2.5 8.0 29.92544 -18.911640 1159321441 Q954
## 2 0 3.0 8.0 26.39530 -14.660804 1159321439 Q953
## 3 0 3.0 8.0 45.87438 15.328226 1159321425 Q805
## 4 0 2.0 7.0 105.38729 21.715416 1159321417 Q881
## 5 0 2.5 7.5 -64.59938 7.182476 1159321411 Q717
## 6 0 5.0 10.0 12.45342 41.903323 1159321407 Q237
## name_ar name_bn name_de name_en name_es
## 1 زيمبابوي জিম্বাবুয়ে Simbabwe Zimbabwe Zimbabue
## 2 زامبيا জাম্বিয়া Sambia Zambia Zambia
## 3 اليمن ইয়েমেন Jemen Yemen Yemen
## 4 فيتنام ভিয়েতনাম Vietnam Vietnam Vietnam
## 5 فنزويلا ভেনেজুয়েলা Venezuela Venezuela Venezuela
## 6 الفاتيكان ভ্যাটিকান সিটি Vatikanstadt Vatican City Ciudad del Vaticano
## name_fa name_fr name_el name_he name_hi name_hu
## 1 زیمبابوه Zimbabwe Ζιμπάμπουε זימבבואה ज़िम्बाब्वे Zimbabwe
## 2 زامبیا Zambie Ζάμπια זמביה ज़ाम्बिया Zambia
## 3 یمن Yémen Υεμένη תימן यमन Jemen
## 4 ویتنام Viêt Nam Βιετνάμ וייטנאם वियतनाम Vietnám
## 5 ونزوئلا Venezuela Βενεζουέλα ונצואלה वेनेज़ुएला Venezuela
## 6 واتیکان Cité du Vatican Βατικανό קריית הוותיקן वैटिकन नगर Vatikán
## name_id name_it name_ja name_ko name_nl name_pl
## 1 Zimbabwe Zimbabwe ジンバブエ 짐바브웨 Zimbabwe Zimbabwe
## 2 Zambia Zambia ザンビア 잠비아 Zambia Zambia
## 3 Yaman Yemen イエメン 예멘 Jemen Jemen
## 4 Vietnam Vietnam ベトナム 베트남 Vietnam Wietnam
## 5 Venezuela Venezuela ベネズエラ 베네수엘라 Venezuela Wenezuela
## 6 Vatikan Città del Vaticano バチカン 바티칸 시국 Vaticaanstad Watykan
## name_pt name_ru name_sv name_tr name_uk name_ur
## 1 Zimbábue Зимбабве Zimbabwe Zimbabve Зімбабве زمبابوے
## 2 Zâmbia Замбия Zambia Zambiya Замбія زیمبیا
## 3 Iémen Йемен Jemen Yemen Ємен یمن
## 4 Vietname Вьетнам Vietnam Vietnam В'єтнам ویتنام
## 5 Venezuela Венесуэла Venezuela Venezuela Венесуела وینیزویلا
## 6 Vaticano Ватикан Vatikanstaten Vatikan Ватикан ویٹیکن سٹی
## name_vi name_zh name_zht fclass_iso tlc_diff fclass_tlc
## 1 Zimbabwe 津巴布韦 辛巴威 Admin-0 country <NA> Admin-0 country
## 2 Zambia 赞比亚 尚比亞 Admin-0 country <NA> Admin-0 country
## 3 Yemen 也门 葉門 Admin-0 country <NA> Admin-0 country
## 4 Việt Nam 越南 越南 Admin-0 country <NA> Admin-0 country
## 5 Venezuela 委内瑞拉 委內瑞拉 Admin-0 country <NA> Admin-0 country
## 6 Thành Vatican 梵蒂冈 梵蒂岡 Admin-0 country <NA> Admin-0 country
## fclass_us fclass_fr fclass_ru fclass_es fclass_cn fclass_tw fclass_in
## 1 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## fclass_np fclass_pk fclass_de fclass_gb fclass_br fclass_il fclass_ps
## 1 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## fclass_sa fclass_eg fclass_ma fclass_pt fclass_ar fclass_jp fclass_ko
## 1 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## fclass_vn fclass_tr fclass_id fclass_pl fclass_gr fclass_it fclass_nl
## 1 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 2 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 3 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 4 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 5 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## 6 <NA> <NA> <NA> <NA> <NA> <NA> <NA>
## fclass_se fclass_bd fclass_ua geometry
## 1 <NA> <NA> <NA> MULTIPOLYGON (((31.28789 -2...
## 2 <NA> <NA> <NA> MULTIPOLYGON (((30.39609 -1...
## 3 <NA> <NA> <NA> MULTIPOLYGON (((53.08564 16...
## 4 <NA> <NA> <NA> MULTIPOLYGON (((104.064 10....
## 5 <NA> <NA> <NA> MULTIPOLYGON (((-60.82119 9...
## 6 <NA> <NA> <NA> MULTIPOLYGON (((12.43916 41...
## [1] 242 169
## [1] "featurecla" "scalerank" "labelrank" "sovereignt" "sov_a3"
## [6] "adm0_dif" "level" "type" "tlc" "admin"
## [11] "adm0_a3" "geou_dif" "geounit" "gu_a3" "su_dif"
## [16] "subunit" "su_a3" "brk_diff" "name" "name_long"
## [21] "brk_a3" "brk_name" "brk_group" "abbrev" "postal"
## [26] "formal_en" "formal_fr" "name_ciawf" "note_adm0" "note_brk"
## [31] "name_sort" "name_alt" "mapcolor7" "mapcolor8" "mapcolor9"
## [36] "mapcolor13" "pop_est" "pop_rank" "pop_year" "gdp_md"
## [41] "gdp_year" "economy" "income_grp" "fips_10" "iso_a2"
## [46] "iso_a2_eh" "iso_a3" "iso_a3_eh" "iso_n3" "iso_n3_eh"
## [51] "un_a3" "wb_a2" "wb_a3" "woe_id" "woe_id_eh"
## [56] "woe_note" "adm0_iso" "adm0_diff" "adm0_tlc" "adm0_a3_us"
## [61] "adm0_a3_fr" "adm0_a3_ru" "adm0_a3_es" "adm0_a3_cn" "adm0_a3_tw"
## [66] "adm0_a3_in" "adm0_a3_np" "adm0_a3_pk" "adm0_a3_de" "adm0_a3_gb"
## [71] "adm0_a3_br" "adm0_a3_il" "adm0_a3_ps" "adm0_a3_sa" "adm0_a3_eg"
## [76] "adm0_a3_ma" "adm0_a3_pt" "adm0_a3_ar" "adm0_a3_jp" "adm0_a3_ko"
## [81] "adm0_a3_vn" "adm0_a3_tr" "adm0_a3_id" "adm0_a3_pl" "adm0_a3_gr"
## [86] "adm0_a3_it" "adm0_a3_nl" "adm0_a3_se" "adm0_a3_bd" "adm0_a3_ua"
## [91] "adm0_a3_un" "adm0_a3_wb" "continent" "region_un" "subregion"
## [96] "region_wb" "name_len" "long_len" "abbrev_len" "tiny"
## [101] "homepart" "min_zoom" "min_label" "max_label" "label_x"
## [106] "label_y" "ne_id" "wikidataid" "name_ar" "name_bn"
## [111] "name_de" "name_en" "name_es" "name_fa" "name_fr"
## [116] "name_el" "name_he" "name_hi" "name_hu" "name_id"
## [121] "name_it" "name_ja" "name_ko" "name_nl" "name_pl"
## [126] "name_pt" "name_ru" "name_sv" "name_tr" "name_uk"
## [131] "name_ur" "name_vi" "name_zh" "name_zht" "fclass_iso"
## [136] "tlc_diff" "fclass_tlc" "fclass_us" "fclass_fr" "fclass_ru"
## [141] "fclass_es" "fclass_cn" "fclass_tw" "fclass_in" "fclass_np"
## [146] "fclass_pk" "fclass_de" "fclass_gb" "fclass_br" "fclass_il"
## [151] "fclass_ps" "fclass_sa" "fclass_eg" "fclass_ma" "fclass_pt"
## [156] "fclass_ar" "fclass_jp" "fclass_ko" "fclass_vn" "fclass_tr"
## [161] "fclass_id" "fclass_pl" "fclass_gr" "fclass_it" "fclass_nl"
## [166] "fclass_se" "fclass_bd" "fclass_ua" "geometry"
# Subset to European countries
europe <- world[world$continent == "Europe", ]
# Display basic information
print(paste("Number of European countries:", nrow(europe)))## [1] "Number of European countries: 50"
## [1] "CRS: WGS 84"
# Create a line connecting cities
city_coords <- st_coordinates(cities_sf)
# Create a LINESTRING
route_matrix <- rbind(
c(-0.1278, 51.5074), # London
c(2.3522, 48.8566), # Paris
c(13.4050, 52.5200) # Berlin
)
route_line <- st_linestring(route_matrix)
route_sf <- st_sf(
route_name = "London-Paris-Berlin",
geometry = st_sfc(route_line, crs = 4326)
)
print(route_sf)## Simple feature collection with 1 feature and 1 field
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -0.1278 ymin: 48.8566 xmax: 13.405 ymax: 52.52
## Geodetic CRS: WGS 84
## route_name geometry
## 1 London-Paris-Berlin LINESTRING (-0.1278 51.5074...
# Create a simple polygon (bounding box around cities)
bbox <- st_bbox(cities_sf)
# Create polygon from bbox
bbox_poly <- st_as_sfc(bbox)
bbox_sf <- st_sf(
name = "Cities Bounding Box",
geometry = bbox_poly
)
print(bbox_sf)## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -3.7038 ymin: 40.4168 xmax: 13.405 ymax: 52.52
## Geodetic CRS: WGS 84
## name geometry
## 1 Cities Bounding Box POLYGON ((-3.7038 40.4168, ...
# Create 100 km buffer around cities
cities_buffer <- st_buffer(cities_sf, dist = 100000) # distance in meters
# Plot
plot(st_geometry(europe), col = "lightgray", main = "Cities with 100km Buffers")
plot(st_geometry(cities_buffer), col = "red", alpha = 0.3, add = TRUE)
plot(st_geometry(cities_sf), col = "blue", pch = 19, add = TRUE)
legend("bottomleft",
legend = c("Countries", "100km Buffer", "Cities"),
col = c("lightgray", "red", "blue"),
pch = c(15, 15, 19))# Find which countries intersect with city buffers
intersecting <- st_intersection(europe, cities_buffer)
# Count intersections per city
intersection_summary <- intersecting %>%
group_by(name.1) %>%
summarise(
countries_within_100km = n(),
.groups = 'drop'
)
print(intersection_summary)## Simple feature collection with 5 features and 2 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4.893563 ymin: 39.50845 xmax: 14.9002 ymax: 53.42907
## Geodetic CRS: WGS 84
## # A tibble: 5 × 3
## name.1 countries_within_100km geometry
## <chr> <int> <POLYGON [°]>
## 1 Berlin 2 ((12.03401 52.18958, 12.02721 52.17413, 12.0415…
## 2 London 1 ((-0.3889138 52.39264, -0.4011157 52.3926, -0.4…
## 3 Madrid 1 ((-4.655691 39.88771, -4.655368 39.88772, -4.65…
## 4 Paris 1 ((2.932405 48.03807, 2.932915 48.04302, 2.96054…
## 5 Rome 2 ((13.54948 41.46206, 13.55526 41.46137, 13.5552…
# Calculate distance matrix between cities (in meters)
dist_matrix <- st_distance(cities_sf)
# Convert to km and create a data frame
dist_km <- units::set_units(dist_matrix, km)
dist_df <- as.data.frame(as.matrix(dist_km))
colnames(dist_df) <- cities_sf$name
rownames(dist_df) <- cities_sf$name
print(round(dist_df, 0))## London Paris Berlin Madrid Rome
## London 0 [km] 344 [km] 932 [km] 1263 [km] 1434 [km]
## Paris 344 [km] 0 [km] 877 [km] 1053 [km] 1105 [km]
## Berlin 932 [km] 877 [km] 0 [km] 1869 [km] 1183 [km]
## Madrid 1263 [km] 1053 [km] 1869 [km] 0 [km] 1364 [km]
## Rome 1434 [km] 1105 [km] 1183 [km] 1364 [km] 0 [km]
# Load libraries
library(sf)
library(dplyr)
library(spData) # spatial example data
# Load country polygons
data("world", package = "spData")
# Load example urban agglomerations (major cities)
data("urban_agglomerations", package = "spData")
# Ensure both layers use the same CRS
urban <- st_transform(urban_agglomerations, st_crs(world))
# Spatial join: city points with country polygons
cities_with_country <- st_join(urban, world, join = st_within)
# Select and rename relevant columns
cities_info <- cities_with_country %>%
select(
city = urban_agglomeration,
country_or_area,
population_millions,
continent, # from world after join
country = name_long
) %>%
st_drop_geometry()
# Print results
print(cities_info)## # A tibble: 540 × 5
## city country_or_area population_millions continent country
## * <chr> <chr> <dbl> <chr> <chr>
## 1 Buenos Aires Argentina 5.17 South America Argentina
## 2 Buenos Aires Argentina 5.91 South America Argentina
## 3 Buenos Aires Argentina 6.76 South America Argentina
## 4 Buenos Aires Argentina 7.55 South America Argentina
## 5 Buenos Aires Argentina 8.42 South America Argentina
## 6 Buenos Aires Argentina 9.14 South America Argentina
## 7 Buenos Aires Argentina 9.92 South America Argentina
## 8 Buenos Aires Argentina 10.5 South America Argentina
## 9 Buenos Aires Argentina 11.1 South America Argentina
## 10 Buenos Aires Argentina 11.8 South America Argentina
## # ℹ 530 more rows
The sp package was the traditional way to handle spatial
data in R. While sf is now preferred, you may encounter
sp objects in older code.
# Note: sp is being phased out, but shown for reference
library(sp)
# Convert sf to sp
cities_sp <- as(cities_sf, "Spatial")
class(cities_sp)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
## [1] "sf" "data.frame"
Raster data represents continuous spatial phenomena as a grid of cells (pixels), each with a value.
The terra package is the modern replacement for the
raster package, offering faster performance and improved
functionality.
# Create a simple raster
r <- rast(nrows = 100, ncols = 100,
xmin = -10, xmax = 10,
ymin = 40, ymax = 60,
crs = "EPSG:4326")
# Assign random values
values(r) <- runif(ncell(r), min = 0, max = 100)
# Display raster properties
print(r)## class : SpatRaster
## size : 100, 100, 1 (nrow, ncol, nlyr)
## resolution : 0.2, 0.2 (x, y)
## extent : -10, 10, 40, 60 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : lyr.1
## min value : 0.009598536
## max value : 99.997116905
## Number of rows: 100
## Number of columns: 100
## Number of cells: 10000
## Resolution: 0.2 0.2
## Extent: -10 10 40 60
## CRS: GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
##
## Summary Statistics:
## lyr.1
## Min. :9.599e-03
## 1st Qu.:2.450e+01
## Median :5.011e+01
## Mean :5.001e+01
## 3rd Qu.:7.527e+01
## Max. :1.000e+02
# Create an elevation-like surface
x_coords <- seq(-10, 10, length.out = 100)
y_coords <- seq(40, 60, length.out = 100)
# Create a matrix with simulated elevation
elevation_matrix <- outer(x_coords, y_coords,
function(x, y) {
100 + 50 * sin(x/2) + 30 * cos(y/3) +
20 * sin(sqrt(x^2 + (y-50)^2)/5)
})
# Create raster from matrix
elevation <- rast(elevation_matrix)
ext(elevation) <- c(-10, 10, 40, 60)
crs(elevation) <- "EPSG:4326"
names(elevation) <- "elevation"
# Plot
plot(elevation,
main = "Simulated Elevation Surface (m)",
col = terrain.colors(100))
contour(elevation, add = TRUE, nlevels = 10)# Create a temperature gradient (decreases with latitude)
temp_values <- matrix(nrow = 100, ncol = 100)
for(i in 1:100) {
for(j in 1:100) {
lat <- 40 + (60-40) * (i-1)/99
temp_values[i, j] <- 25 - (lat - 40) * 0.5 + rnorm(1, 0, 1)
}
}
temperature <- rast(temp_values)
ext(temperature) <- c(-10, 10, 40, 60)
crs(temperature) <- "EPSG:4326"
names(temperature) <- "temperature"
# Plot
plot(temperature,
main = "Simulated Temperature in degree Celsius",
col = rev(heat.colors(100)))# Create multiple rasters
r1 <- elevation
r2 <- temperature
# Basic arithmetic
r_sum <- r1 + r2
r_diff <- r1 - r2
r_product <- r1 * r2
# Plot comparison
par(mfrow = c(2, 2))
plot(r1, main = "Elevation", col = terrain.colors(50))
plot(r2, main = "Temperature", col = heat.colors(50))
plot(r_sum, main = "Sum", col = viridis(50))
plot(r_diff, main = "Difference", col = viridis(50))# Classify elevation into categories
# Low: < 110, Medium: 110-130, High: > 130
rcl_matrix <- matrix(c(
-Inf, 110, 1,
110, 130, 2,
130, Inf, 3
), ncol = 3, byrow = TRUE)
elevation_class <- classify(elevation, rcl_matrix)
# Create labels
levels(elevation_class) <- data.frame(
value = 1:3,
category = c("Low", "Medium", "High")
)
# Plot
plot(elevation_class,
main = "Elevation Classification",
col = c("green", "yellow", "brown"))# Calculate mean in 3x3 moving window
elevation_smooth <- focal(elevation, w = 3, fun = mean)
# Compare original and smoothed
par(mfrow = c(1, 2))
plot(elevation, main = "Original Elevation", col = terrain.colors(100))
plot(elevation_smooth, main = "Smoothed Elevation (3x3)", col = terrain.colors(100))# Define a smaller extent
crop_ext <- ext(-5, 5, 45, 55)
# Crop raster
elevation_crop <- crop(elevation, crop_ext)
# Create a circular mask
center_x <- 0
center_y <- 50
radius <- 5
# Create mask polygon
circle <- st_buffer(st_point(c(center_x, center_y)), dist = radius)
circle_sf <- st_sf(geometry = st_sfc(circle, crs = 4326))
# Mask raster
elevation_mask <- mask(elevation_crop, vect(circle_sf))
# Plot
par(mfrow = c(1, 2))
plot(elevation_crop, main = "Cropped", col = terrain.colors(100))
plot(elevation_mask, main = "Masked (Circular)", col = terrain.colors(100))# Extract raster values at city locations
cities_elevation <- extract(elevation, vect(cities_sf))
cities_temperature <- extract(temperature, vect(cities_sf))
# Combine with city data
cities_with_raster <- cities_sf %>%
mutate(
elevation_m = cities_elevation$elevation,
temperature_c = cities_temperature$temperature
) %>%
st_drop_geometry()
print(cities_with_raster)## name population elevation_m temperature_c
## 1 London 8982000 53.56876 21.18578
## 2 Paris 2161000 141.07056 20.14503
## 3 Berlin 3645000 NA NA
## 4 Madrid 3223000 39.06862 15.33385
## 5 Rome 2873000 NA NA
# Create a multi-layer raster (stack)
environmental_stack <- c(elevation, temperature)
names(environmental_stack) <- c("elevation", "temperature")
print(environmental_stack)## class : SpatRaster
## size : 100, 100, 2 (nrow, ncol, nlyr)
## resolution : 0.2, 0.2 (x, y)
## extent : -10, 10, 40, 60 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## names : elevation, temperature
## min values : 34.38923, 11.76533
## max values : 199.86452, 27.60590
# Plot all layers
plot(environmental_stack,
col = viridis(100),
main = c("Elevation (m)", "Temperature in degree Celsius"))Understanding coordinate reference systems (CRS) is crucial for accurate geospatial analysis.
A Coordinate Reference System defines how spatial data relates to locations on Earth.
## Cities CRS:
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## MEMBER["World Geodetic System 1984 (G2296)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
##
## World CRS:
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
##
## Same CRS? TRUE
EPSG codes are standardized numerical identifiers for CRS.
# Create a reference table
epsg_reference <- data.frame(
EPSG = c(4326, 3857, 27700, 32630, 2154, 3035),
Name = c("WGS 84", "Web Mercator", "British National Grid",
"WGS 84 / UTM zone 30N", "RGF93 / Lambert-93",
"ETRS89 / LAEA Europe"),
Type = c("Geographic", "Projected", "Projected",
"Projected", "Projected", "Projected"),
Use_Case = c("GPS, Global", "Web maps", "UK",
"Western Europe", "France", "Europe-wide")
)
print(epsg_reference)## EPSG Name Type Use_Case
## 1 4326 WGS 84 Geographic GPS, Global
## 2 3857 Web Mercator Projected Web maps
## 3 27700 British National Grid Projected UK
## 4 32630 WGS 84 / UTM zone 30N Projected Western Europe
## 5 2154 RGF93 / Lambert-93 Projected France
## 6 3035 ETRS89 / LAEA Europe Projected Europe-wide
## EPSG 4326 Details:
## Name: WGS 84
## Datum: WGS84
## Units: degree
## Is geographic? TRUE
PROJ4 strings are another way to define CRS, though EPSG codes are preferred.
# Get PROJ4 string
proj4_wgs84 <- st_crs(4326)$proj4string
cat("EPSG:4326 PROJ4 string:\n", proj4_wgs84, "\n\n")## EPSG:4326 PROJ4 string:
## +proj=longlat +datum=WGS84 +no_defs
# Create custom projection
custom_crs <- st_crs("+proj=aeqd +lat_0=51.5 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m")
cat("Custom azimuthal equidistant projection:\n")## Custom azimuthal equidistant projection:
## Coordinate Reference System:
## User input: +proj=aeqd +lat_0=51.5 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m
## wkt:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["unknown",
## METHOD["Azimuthal Equidistant",
## ID["EPSG",1125]],
## PARAMETER["Latitude of natural origin",51.5,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]
# Transform cities to different projections
cities_webmercator <- st_transform(cities_sf, crs = 3857) # Web Mercator
cities_utm <- st_transform(cities_sf, crs = 32630) # UTM Zone 30N
# Compare coordinates
comparison <- data.frame(
City = cities_sf$name,
WGS84_lon = st_coordinates(cities_sf)[,1],
WGS84_lat = st_coordinates(cities_sf)[,2],
WebMerc_x = st_coordinates(cities_webmercator)[,1],
WebMerc_y = st_coordinates(cities_webmercator)[,2],
UTM_x = st_coordinates(cities_utm)[,1],
UTM_y = st_coordinates(cities_utm)[,2]
)
print(comparison)## City WGS84_lon WGS84_lat WebMerc_x WebMerc_y UTM_x UTM_y
## 1 London -0.1278 51.5074 -14226.63 6711542 699316.2 5710164
## 2 Paris 2.3522 48.8566 261845.71 6250564 892519.3 5425340
## 3 Berlin 13.4050 52.5200 1492237.77 6894700 1608986.4 5946385
## 4 Madrid -3.7038 40.4168 -412305.13 4926697 440290.5 4474257
## 5 Rome 12.4964 41.9028 1391092.88 5146430 1786888.1 4756736
# Transform elevation raster to UTM
elevation_utm <- project(elevation, "EPSG:32630", method = "bilinear")
# Compare
cat("Original CRS:", crs(elevation, describe = TRUE)$name, "\n")## Original CRS: WGS 84
## Transformed CRS: WGS 84 / UTM zone 30N
# Plot comparison
par(mfrow = c(1, 2))
plot(elevation, main = "WGS84 (Geographic)", col = terrain.colors(100))
plot(elevation_utm, main = "UTM Zone 30N (Projected)", col = terrain.colors(100))CRS choice affects area and distance calculations significantly.
# Calculate area of European countries in different CRS
europe_wgs84 <- europe
europe_laea <- st_transform(europe, crs = 3035) # LAEA Europe
# Calculate areas
area_wgs84 <- st_area(europe_wgs84)
area_laea <- st_area(europe_laea)
# Compare for a few countries
comparison_df <- data.frame(
Country = europe$admin[1:5],
Area_WGS84_km2 = as.numeric(area_wgs84[1:5]) / 1e6,
Area_LAEA_km2 = as.numeric(area_laea[1:5]) / 1e6
) %>%
mutate(Difference_pct = abs(Area_WGS84_km2 - Area_LAEA_km2) / Area_LAEA_km2 * 100)
print(comparison_df)## Country Area_WGS84_km2 Area_LAEA_km2 Difference_pct
## 1 Vatican 7.064538e-01 7.075141e-01 0.1498622
## 2 Jersey 1.263326e+02 1.267385e+02 0.3203309
## 3 Guernsey 4.805537e+01 4.821272e+01 0.3263626
## 4 Isle of Man 5.573146e+02 5.597506e+02 0.4351966
## 5 United Kingdom 2.396136e+05 2.406448e+05 0.4285199
cat("\nNote: WGS84 (geographic) gives incorrect areas. Always use projected CRS for area calculations!\n")##
## Note: WGS84 (geographic) gives incorrect areas. Always use projected CRS for area calculations!
Creating publication-quality static maps is essential for reports and papers.
# Simple map of Europe
ggplot(data = europe) +
geom_sf(fill = "lightblue", color = "white", size = 0.2) +
coord_sf(xlim = c(-25, 45), ylim = c(35, 72)) +
theme_minimal() +
labs(title = "Map of Europe",
subtitle = "Using ggplot2 and geom_sf")# Map with cities
ggplot() +
geom_sf(data = europe, fill = "gray90", color = "white", size = 0.3) +
geom_sf(data = cities_sf, aes(size = population),
color = "red", alpha = 0.6) +
geom_sf_text(data = cities_sf, aes(label = name),
nudge_y = 1, size = 3, fontface = "bold") +
scale_size_continuous(name = "Population",
labels = scales::comma,
range = c(2, 10)) +
coord_sf(xlim = c(-10, 20), ylim = c(40, 60)) +
theme_minimal() +
labs(title = "Major European Cities",
subtitle = "Population size indicated by point size",
x = "Longitude", y = "Latitude")# Map population density
europe_density <- europe %>%
mutate(pop_density = pop_est / (as.numeric(st_area(.)) / 1e6)) # per km²
ggplot(data = europe_density) +
geom_sf(aes(fill = pop_density), color = "white", size = 0.2) +
scale_fill_viridis_c(name = "Population Density\n(per sq KM)",
trans = "log10",
labels = scales::comma) +
coord_sf(xlim = c(-25, 45), ylim = c(35, 72)) +
theme_minimal() +
labs(title = "Population Density in Europe",
subtitle = "Log scale transformation applied")# Complex map with multiple layers
ggplot() +
geom_sf(data = europe, fill = "gray95", color = "gray70", size = 0.3) +
geom_sf(data = cities_buffer, fill = "orange", alpha = 0.1, color = "orange") +
geom_sf(data = route_sf, color = "blue", size = 1.5, linetype = "dashed") +
geom_sf(data = cities_sf, aes(size = population),
color = "red", alpha = 0.8) +
geom_sf_text(data = cities_sf, aes(label = name),
nudge_y = 1.5, size = 3.5, fontface = "bold") +
scale_size_continuous(name = "Population",
labels = scales::comma,
range = c(3, 12)) +
coord_sf(xlim = c(-12, 18), ylim = c(40, 58)) +
theme_minimal() +
theme(legend.position = "bottom") +
labs(title = "European Cities with Travel Route and Buffers",
subtitle = "100km buffers shown in orange",
x = "Longitude", y = "Latitude")# Create categories for population
europe_categories <- europe %>%
dplyr::filter(continent == "Europe") %>%
mutate(pop_category = cut(pop_est,
breaks = c(0, 5e6, 20e6, 50e6, Inf),
labels = c("< 5M", "5M - 20M", "20M - 50M", "> 50M")))
ggplot(data = europe_categories) +
geom_sf(aes(fill = pop_category), color = "white", size = 0.2) +
scale_fill_viridis_d(name = "Population") +
facet_wrap(~ pop_category) +
theme_minimal() +
theme(axis.text = element_blank(),
axis.ticks = element_blank()) +
labs(title = "European Countries by Population Category")tmap provides a grammar of graphics specifically
designed for maps.
# Simple map
tm_shape(europe) +
tm_borders(col = "gray50") +
tm_fill(col = "lightblue") +
tm_layout(title = "Europe - Basic tmap")# Population map
tm_shape(europe) +
tm_polygons("pop_est",
title = "Population",
palette = "YlOrRd",
style = "jenks",
n = 5) +
tm_layout(title = "European Population",
legend.outside = TRUE,
legend.outside.position = "right",
frame = FALSE)# Complex map with multiple layers
tm_shape(europe) +
tm_polygons("pop_est",
title = "Population",
palette = "Blues",
style = "quantile",
border.col = "white",
border.alpha = 0.5) +
tm_shape(cities_sf) +
tm_symbols(size = "population",
col = "red",
alpha = 0.7,
title.size = "City Population",
scale = 1.5) +
tm_text("name",
size = 0.7,
auto.placement = TRUE,
fontface = "bold") +
tm_layout(title = "Europe: Country and City Populations",
legend.outside = TRUE,
legend.outside.position = "right",
bg.color = "lightcyan",
frame = TRUE) +
tm_compass(position = c("left", "top")) +
tm_scale_bar(position = c("left", "bottom"))# Select a few countries for detailed view
selected_countries <- c("France", "Germany", "Spain", "Italy")
countries_subset <- europe %>%
dplyr::filter(admin %in% selected_countries)
# Faceted map
tm_shape(countries_subset) +
tm_polygons("pop_est",
title = "Population",
palette = "Greens") +
tm_facets(by = "admin",
free.coords = TRUE,
ncol = 2) +
tm_layout(legend.show = FALSE,
panel.labels = selected_countries)# Map elevation raster
tm_shape(elevation) +
tm_raster(title = "Elevation (m)",
palette = terrain.colors(100),
style = "cont") +
tm_shape(cities_sf) +
tm_symbols(size = 0.5, col = "red") +
tm_text("name", size = 0.7, auto.placement = TRUE) +
tm_layout(title = "Elevation with Cities",
legend.outside = TRUE,
bg.color = "lightblue") +
tm_compass(position = c("right", "top")) +
tm_scale_bar(position = c("left", "bottom"))Interactive maps allow users to zoom, pan, and click on features for more information.
# Create custom popups
popup_content <- paste0(
"<strong>", cities_sf$name, "</strong><br/>",
"Population: ", format(cities_sf$population, big.mark = ","), "<br/>",
"Coordinates: ", round(st_coordinates(cities_sf)[,2], 2), "°N, ",
round(st_coordinates(cities_sf)[,1], 2), "°E"
)
# Ensure UTF-8
cities_sf$name <- iconv(cities_sf$name, from = "", to = "UTF-8", sub = "")
popup_content <- paste0("<b>", cities_sf$name, "</b><br>Population: ", cities_sf$population)
popup_content <- iconv(popup_content, from = "", to = "UTF-8", sub = "")
# Map with custom popups and circle markers
leaflet(cities_sf) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addCircleMarkers(
radius = ~sqrt(population) / 100,
color = "red",
fillColor = "orange",
fillOpacity = 0.6,
popup = popup_content,
label = ~name
) %>%
addLegend(
position = "bottomright",
colors = "orange",
labels = "European Cities",
title = "Legend"
)# Create color palette
pal <- colorNumeric(palette = "YlOrRd", domain = europe$pop_est)
# Choropleth map
leaflet(europe) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
fillColor = ~pal(pop_est),
fillOpacity = 0.7,
color = "white",
weight = 1,
popup = ~paste0("<strong>", admin, "</strong><br/>",
"Population: ", format(pop_est, big.mark = ",")),
highlight = highlightOptions(
weight = 3,
color = "red",
fillOpacity = 0.9,
bringToFront = TRUE
)
) %>%
addLegend(
pal = pal,
values = ~pop_est,
title = "Population",
position = "bottomright",
labFormat = labelFormat(big.mark = ",")
)# Color palette for raster
raster_pal <- colorNumeric(terrain.colors(100),
values(elevation),
na.color = "transparent")
# Complex map with multiple layers
leaflet() %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
addProviderTiles(providers$OpenStreetMap, group = "Streets") %>%
addRasterImage(raster(elevation),
colors = raster_pal,
opacity = 0.6,
group = "Elevation") %>%
addPolygons(data = europe,
fillColor = "transparent",
color = "blue",
weight = 2,
group = "Countries") %>%
addCircleMarkers(data = cities_sf,
radius = ~sqrt(population) / 100,
color = "red",
fillColor = "yellow",
fillOpacity = 0.8,
popup = ~paste0("<strong>", name, "</strong><br/>",
"Pop: ", format(population, big.mark = ",")),
group = "Cities") %>%
addPolylines(data = route_sf,
color = "purple",
weight = 3,
opacity = 0.7,
group = "Route") %>%
addLayersControl(
baseGroups = c("Satellite", "Streets"),
overlayGroups = c("Elevation", "Countries", "Cities", "Route"),
options = layersControlOptions(collapsed = FALSE)
) %>%
addLegend(
pal = raster_pal,
values = values(elevation),
title = "Elevation (m)",
position = "bottomright"
)mapview provides a simple interface for quick
interactive visualization.
# More control over appearance
mapview(cities_sf,
zcol = "population",
col.regions = c("yellow", "orange", "red"),
at = c(0, 3000000, 6000000, 9000000),
cex = 6,
alpha = 0.8,
legend = TRUE,
layer.name = "City Population",
map.types = c("OpenStreetMap", "Esri.WorldImagery"),
popup = leafpop::popupTable(
st_drop_geometry(cities_sf),
zcol = c("name", "population")
))# Calculate which areas are within 200km of major cities
library(sf)
library(dplyr)
library(ggplot2)
library(viridis)
# --------------------------
# Step 0: Prepare data
# --------------------------
# cities_sf = your 5 major European cities
cities_sf <- st_sf(
name = c("London", "Paris", "Berlin", "Madrid", "Rome"),
population = c(8982000, 2161000, 3645000, 3223000, 2873000),
geometry = st_sfc(
st_point(c(-0.1278, 51.5074)),
st_point(c(2.3522, 48.8566)),
st_point(c(13.405, 52.52)),
st_point(c(-3.7038, 40.4168)),
st_point(c(12.4964, 41.9028))
),
crs = 4326
)
# europe = your European country polygons
# europe <- st_read("path_to_europe_shapefile.shp")
# --------------------------
# Step 1: Create 200 km buffers
# --------------------------
cities_200km <- st_transform(cities_sf, 3035) %>%
st_buffer(200000)
cities_200km_wgs84 <- st_transform(cities_200km, 4326)
# --------------------------
# Step 2: Calculate coverage
# --------------------------
coverage <- st_union(cities_200km_wgs84)
# --------------------------
# Step 3: Find countries with >50% coverage
# --------------------------
europe_3035 <- st_transform(europe, 3035)
cities_200km_3035 <- st_transform(cities_200km_wgs84, 3035)
coverage_analysis <- europe_3035 %>%
rowwise() %>%
mutate(
country_area = as.numeric(st_area(geometry)),
covered_area = {
intr <- st_intersection(geometry, st_union(cities_200km_3035))
intr <- st_collection_extract(intr, "POLYGON") # only polygons
if (length(intr) == 0) 0 else as.numeric(st_area(st_union(intr)))
},
coverage_pct = (covered_area / country_area) * 100
) %>%
ungroup() %>%
st_transform(4326) %>%
arrange(desc(coverage_pct))
# --------------------------
# Step 4: Visualize results
# --------------------------
ggplot() +
geom_sf(data = coverage_analysis, aes(fill = coverage_pct),
color = "white", size = 0.2) +
geom_sf(data = cities_200km_wgs84, fill = "transparent", color = "red",
linetype = "dashed", alpha = 0.3) +
geom_sf(data = cities_sf, color = "red", size = 3) +
scale_fill_viridis_c(name = "Coverage %", option = "plasma") +
theme_minimal() +
labs(title = "Accessibility: Areas within 200km of Major Cities",
subtitle = "Five major European cities analyzed")# --------------------------
# Step 5: Print top results
# --------------------------
cat("\nTop 10 Countries by Coverage:\n")##
## Top 10 Countries by Coverage:
## # A tibble: 10 × 2
## admin coverage_pct
## <chr> <dbl>
## 1 Vatican 100.
## 2 United Kingdom 29.8
## 3 Germany 25.3
## 4 Spain 24.8
## 5 Italy 20.8
## 6 France 19.1
## 7 Poland 10.1
## 8 Belgium 2.33
## 9 Czechia 0.678
## 10 Jersey 0
# Identify high-elevation, low-temperature areas
# Step 1: Define thresholds
high_elevation <- elevation > 130
low_temperature <- temperature < 15
# Step 2: Combine conditions
suitable_areas <- high_elevation & low_temperature
# Step 3: Vectorize suitable areas
suitable_polygons <- as.polygons(suitable_areas) %>%
st_as_sf() %>%
dplyr::filter(elevation == 1) # Keep only TRUE values
# Step 4: Calculate area
suitable_polygons_3035 <- st_transform(suitable_polygons, crs = 3035)
total_area <- sum(as.numeric(st_area(suitable_polygons_3035))) / 1e6 # km²
# Step 5: Visualize
ggplot() +
geom_sf(data = suitable_polygons, fill = "purple", alpha = 0.5) +
geom_sf(data = cities_sf, color = "red", size = 3) +
theme_minimal() +
labs(title = "Suitable Areas: High Elevation (>130m) & Low Temperature (<15 degree Celsius)",
subtitle = paste0("Total suitable area: ", round(total_area, 0), " sq KM"))##
## Suitable area statistics:
## Total area: 5678.26 km²
## Number of polygons: 1
# 1. Use appropriate CRS for analysis
# - Geographic (EPSG:4326) for global visualization
# - Projected CRS for measurements and analysis
# 2. Simplify geometries when appropriate
simplified <- st_simplify(europe, dTolerance = 1000)
# 3. Use spatial indexing for large datasets
# st_make_valid() # Fix invalid geometries
# st_is_valid() # Check validity
# 4. For rasters, use appropriate resolution
# Aggregate to coarser resolution for faster processing
elevation_coarse <- aggregate(elevation, fact = 5, fun = mean)
# 5. Cache intermediate results
# saveRDS(large_sf_object, "cache/processed_data.rds")
# large_sf_object <- readRDS("cache/processed_data.rds")# 1. NEVER calculate areas in geographic CRS (EPSG:4326)
# WRONG:
area_wrong <- st_area(europe) # If europe is in EPSG:4326
# RIGHT:
europe_projected <- st_transform(europe, crs = 3035)
area_correct <- st_area(europe_projected)
# 2. Always check CRS compatibility
# st_crs(layer1) == st_crs(layer2)
# 3. Handle invalid geometries
# layer_valid <- st_make_valid(layer_with_issues)
# 4. Be careful with st_intersection on large datasets
# Use st_crop first to reduce computation
# cropped <- st_crop(large_layer, small_extent)
# result <- st_intersection(cropped, other_layer)
# 5. Remember raster vs vector differences
# Rasters: continuous data, fixed resolution
# Vectors: discrete features, variable precision# Export vector data
st_write(cities_sf, "output/cities.shp") # Shapefile
st_write(cities_sf, "output/cities.geojson") # GeoJSON
st_write(cities_sf, "output/cities.gpkg") # GeoPackage (recommended)
# Export raster data
writeRaster(elevation, "output/elevation.tif") # GeoTIFF
# Export maps
ggsave("output/map.png", width = 10, height = 8, dpi = 300)
tmap_save(tm_map_object, "output/tmap.png", dpi = 300)| Package | Purpose | Key Functions |
|---|---|---|
sf |
Vector data | st_read(), st_transform(),
st_buffer(), st_intersection() |
terra |
Raster data | rast(), project(), crop(),
mask(), extract() |
ggplot2 |
Static maps | geom_sf(), coord_sf() |
tmap |
Thematic maps | tm_shape(), tm_polygons(),
tm_raster() |
leaflet |
Interactive maps | leaflet(), addTiles(),
addPolygons() |
mapview |
Quick viewing | mapview() |
This material is part of the training program by The National Centre for Research Methods © NCRM authored by Dr Somnath Chaudhuri (University of Southampton). Content is under a CC BY‑style permissive license and can be freely used for educational purposes with proper attribution.